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Abstract 

Bound-state solutions are obtained numerically in the instantaneous ap- 
proximation for a spin-0 and spin-1/2 constituent that interact via minimal 
electrodynamics. To solve the integral equations in momentum space, a 
method is developed for integrating over the logarithmic singularity in ker- 
nels, making it possible to use basis functions that essentially automatically 
satisfy the boundary conditions. For bound-state solutions that decrease 
rapidly at small and large values of momentum, accurate solutions are ob- 
tained with significantly fewer basis functions when the solution is expanded 
in terms of these more general basis functions. The presence of a deriva- 
tive coupling in single-photon exchange complicates the construction of the 
Bethe-Salpeter equation in the instantaneous approximation and, in the non- 
relativistic limit, gives rise to an additional electrostatic potential term that 
is second order in the coupling constant and decreases as the square of the 
distance between constituents. 



PACS numbers: 02.60.-x, 03.65.Pm, 11. 10. St 



1 Introduction 



There has been a growing interest in solving relativistic, bound-state 
equations both because important bound-state systems are relativistic and 
because the development of high-speed computers makes it possible to solve 
such equations. The Bcthc-Salpctcr equation [1], which is based on field the- 
ory, is covariant and reduces to the Schrodinger equation in the nonrelativis- 
tic limit, is the appropriate equation to use in describing relativistic bound 
states. Unfortunately, even numerically the two-body bound-state equation 
is exceedingly difficult to solve [2]. For this reason various approximations 
such as the Blankenbecler-Sugar approximation [3] or the instantaneous ap- 
proximation[l,4] are often made that reduce the covariant equation in four di- 
mensions to an approximately-covariant equation in three dimensions. In this 
article attention is restricted to the instantaneous approximation although 
the general methods developed here can be used in the implementation of 
other approximation schemes that reduce the Bethe-Salpeter equation to a 
three-dimensional equation. 

The instantaneous approximation, which is the approximation that the 
binding quanta travel instantaneously between the bound constituents, was 
first introduced in the original article by Bethe and Salpetcr [1] and was 
used in their calculation to demonstrate that the Schrodinger equation is 
the nonrelativistic limit of the Bethe-Salpeter equation. More in the spirit 
of this work, Salpeter [4] made the instantaneous approximation to reduce 
the Bethe-Salpeter equation to a three-dimensional equation and calculated 
corrections to the fine structure of hydrogen-like atoms. The presence of a 
derivative coupling in single-photon exchange complicates the construction 
of the Bethe-Salpeter equation in the instantaneous approximation and, in 
addition to terms that are first order in the coupling constant, gives rise to two 
terms that are second order. The "seagull" interaction yields the same two 
second-order interactions, but with different strengths. In the nonrelativistic 
limit, the second-order interaction becomes an electrostatic potential term 
that decreases as the square of the distance between the constituents. 

A major difficulty in solving equations numerically when the instanta- 
neous approximation is made arises because a logarithmic singularity occurs 
in the kernel of the integral equations. Gammel and Menzel [5] overcome the 
problem by using a special weighting scheme in the neighborhood of the sin- 
gularity, and Eyre and Vary [6] introduce a numerical cutoff and then correct 
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for the effects of the cutoff using perturbation theory. Later Spence and Vary 
[7] use B-sphnes [8] as basis functions and perform all integrals analytically. 
Since the B-splines are polynomials, their method is restricted to polynomial 
basis functions. In this article a method is used to integrate over the loga- 
rithmic singularity that allows the use of more general basis functions that 
essentially automatically satisfy the boundary conditions and are not neces- 
sarily polynomials. This new method, which is very simple conceptually, is 
a generalization of the method introduced in Ref. 7. For bound states that 
decrease rapidly at small and large values of momentum, accurate solutions 
are obtained with significantly fewer basis functions when these more general 
basis functions are used. 

Two physical problems that are of immediate interest are constituent 
models of quarks and leptons [9,10] and constituent-quark models of mesons 
[11]. Equations that account for some relativistic effects have had success in 
describing the properties of both light and heavy mesons [11]. The fact that 
there are three families of leptons, much as there are families of elements 
or families of hadrons, suggests that the leptons might be composite. The 
discovery of neutrino oscillations [12] raises the possibility that neutrinos 
might also be composite [13]. If the electron, muon and tau are bound 
states of a single system, the system is necessarily relativistic: The mass of 
the tau must, of course, be less than the sum of the masses of the bound, 
constituent particles. Since the ratio of the electron's mass to that of the 
tau's equals 1/3536, the mass of the electron is less than 1/3536 of the sum 
of the constituent masses, indicating highly relativistic binding. 

The Bethe-Salpeter equation discussed here has been solved exactly in the 
strong binding (zero energy) limit when the ladder approximation is made 
[10]. In the instantaneous approximation the numerical solutions obtained 
here are not in good agreement with the exact zero-energy solutions so, not 
surprisingly, the instantaneous approximation is not satisfactory for very 
strongly bound states. 

To estimate the accuracy of each solution, in the physical region the left- 
and right-hand sides of the equation are calculated midway between each 
knot, and a rehabihty coefficient R [14], which is a statistical measure of how 
accurately the left- and right-hand sides agree at the selected points, is calcu- 
lated. Examining points where the left- and right-hand sides of the equation 
agree least well reveals possible problems with solutions and suggests possible 
remedies. 
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2 Derivation and Separation of the Bethe— 
Salpeter Equation in the Instantaneous 
Approximation 

When a spin-0 field (f){x), which represents a quanta with charge Q and 
mass M, interacts via minimal electrodynamics with a spin- 1/2 field ^(x), 
which represents a quanta with charge q and mass m, the renormalizable 
Lagrangian is [15] 

L =: [{id'^ - QA'^Mi-id^ - QA^)<P+] - M^+cl> 

+ ^'7^,(^(9^ - qA^)^> - m^* - ^F^^F^"" : (2.1) 
where F^,, = - 5^^^. 

The two-particle, Bethe-Salpeter wave function is defined by 

Xk{xuX2) =< ^\T{^{x^)^{x2))\K > . (2.2) 

In (2.2) the symbol T represents time ordering and the letter K labels the 
four-momentum of the bound state. The center-of-mass coordinates are 
defined by 

= ixt + (1 - e)x^, (2.3) 
and the relative coordinates x'^ by 

x^^xl- x^. (2.4) 

When the parameter ^ is given by ^ = m/{m -\- M) , the usual nonrela- 
tivistic definition of center-of-mass coordinates results. As will be seen, the 
parameter C, drops out of the Bethe-Salpeter equation when the instanta- 
neous approximation is made so there is no need to make a specific choice. 
The dependence of xk{xi, X2) on the center-of-mass coordinates factors with 
the result that Xk{xi,X2) can be rewritten as 

Xk{xi, X2) = (27r)-=^/2e-^^^''Xi.(x). (2.5) 
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Denoting the Fourier transform of Xk{x) by XKip), the Bethe-Salpeter 
equation is 



(p^^ + Cii-'^TM - ^){b^ - (1 - O^lb/. - (1 - OK,] - m'}xk{p) 

iqQ d^q 



(27r)4 J-oo (p - g)2 + ie 
?7r)8 J- 



A{qQy (fq 2m - ^''7^ 

^ (27r)8 i-oo g2 _ + (p - g + ^ii:)2 + ie ^ 



/ 



-Xi.(^). (2.6) 



(g - A; - ^Ky + ie' 

A charged, spin-0 boson interacts electromagnetically through two funda- 
mentally different processes: single-photon exchange and the "seagull" in- 
teraction. In the above equation the terms proportional to qQ and {qQY 
arise, respectively, from these two interactions. Although the "seagull" in- 
teraction is second order in the coupling constant, it has been included in 
the Bethe-Salpeter equation to determine its effect on solutions. The spirit 
of this calculation, then, is similar to others dealing with the Bethe-Salpeter 
equation whereby the ladder approximation is made, but some solutions with 
large coupling constants are studied. The Feynman diagrams for these two 
interactions are shown in Fig. 2.1 

The instantaneous approximation [1,4] is made by making the replace- 
ment 

11 1 

(2.7) 



1 1 1 



+ ie /cg ~ + -|- ie 

in each of the three photon propagators in (2.6). Defining 

/oo 
dkoXKik), (2.8a) 
-00 

and 

/oo 
dkokoxxik), (2.86) 
-00 

the integral over go in the single-photon-exchange term in (2.6) and the in- 
tegral over ko in the "seagull" term can be carried out immediately. In the 
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Figure 2.1: Feynman diagrams for (a) single-photon exchange and (b) the 
"seagull" interaction. The sohd, dashed and wavy hues represent a spin- 1/2 
fermion, a spin-0 boson and a photon, respectively. 

"seagull" term it is then possible to integrate over qq. Going to the center- 
of-mass frame where = {E, 0), (2.6) becomes 

- pY + ^E^o - m]{\p' - (1 - 0^]' - pV - M'}xEip) 

iqQ d?q 



(27r)4 y_oo (p - q; 
iqQ d?q 



(27r)4 y-oo (p - q)2 
2i{qQf r°° d^q 2m + 



r ur^^^y - ^'(P' + ?0 - 2(1 - 0^/]*^(q) 

J—oo IP — QJ 

7VE(q) 

/oo d'^ fc 
-oo(^^ ^^-^^ 

In the above equation Xe{p) is the value of Xi^(p) in the center-of-mass frame, 
etc., and a;^(g) = (m^ + g^)^/^. 

Solving (2.9) for Xsip), integrating over p° on both sides of the equation 
and using (2.8) yields 

u;m{p)E^e{p) = ^^[u;M{p)+u;m{p)]b'Yp' + ^'rn]^E{p) 



5 



i-oo p - q 2 a;„, « 



167r3 i-oo (p - q)2 a;^(p) 

>E(q) 



167r''^ J-oo (p — q)2 



(27r)6 i-oo (p - q)2 i-oo (q - k)^ ^ " ^ 

Eq. (2.10) cannot readily be solved because of the presence of the two func- 
tions '^E and that are related as indicated in (2.8). The function 0^ is 
present, of course, as a consequence of the derivative coupling. 

It is possible, however, to express 0^; in terms of e as follows: By 
dividing both sides of (2.9) by — (1 — 0-^]^ ~ PV — M"^} and then 
integrating over a second equation is obtained that involves both and 

(t>E'- 

qQ f°° (fq 



+ 



r [(i-0^ + 7V(P^ + gO]Mq) 

j-oo \p — q; 

.(q) 



167r3 

qQ d^q 



IGtt^ i-oo (p — q)2 

{qQY f°° (fq 2m7° + 7V?* f°° d^k 



a^q 2mr + 7"7 g r m (^.\ 



(27r)6 

Subtracting (2.11) from (2.10) and solving for 4>e{p) yields 

0e(p) = (i-O£^*b(p) 
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1^ (/,V + /„)[.„(rt*,(p) + ^/~^*,«,)|. (2.12) 



Using (2.12) to express 0£:(p) in terms of (2.10) becomes the Bethe- 
Salpeter equation in the instantaneous approximation: 



qQ 



167r3 



(p - q)' 



,'^m(p) 



7%V + 



,'^m(?) 



+ 



7°m ^ ^i5(q) 



.^m(p) '^m(g). 



4(27r) 



00 (p-qj 



(q - k): 



r*E(k) 



(9Q) /""^ 2m7" + 7^7* (Fk 



(27r)6 J-00 (p-q)2 a;„(g) 



(q-k)^ 



,*E(k) (2.13) 



In (2.13) the final two terms, which are proportional to {qQY, arise from the 
derivative coupling in single photon exchange and the "seagull" interaction, 
respectively. 

Eq. (2.13) is much easier to solve numerically than the Bethe-Salpeter 
equation, primarily because it is much easier to obtain solutions with real 
energy eigenvalues. Specifically, equations of the form (2.13) can be solved 
numerically by converting them to matrix eigenvalue equations. When each 
side is multiphed by ^|;(p) and integrated over dPp, excluding the eigenvalue 

the quantity on the left-hand side is Hermitian and positive definite and 
the quantity on the right-hand side is Hermitian. As a consequence the 
energy eigenvalue must be real [16]. The Hermiticity results from the fact 
that the momenta p and q appear symmetrically on the right-hand side of 
(2.13). In the very special cases where the Bethe-Salpeter equation possesses 
the Hermiticity properties of (2.13), the equation is relatively easy to solve 
numerically [2,17] in spite of the fact that after separation of the two angular 
variables, it is still an integral or partial differential equation in two variables. 



7 



Solutions to (2.13) are of the form 



(2.14) 



where the (t)^^\6^Lp) are the same functions [15] that represent the angular 
dependence of the bound-state solutions to the Dirac equation when the po- 
tential is spherically symmetric. After the angular integration is performed 
using Hecke's theorem [18,19,20,10], the angular variables separate. Multi- 
plying the resulting upper and lower equations by p and — p, respectively, 
yields 



Eujm{p) 



[uJM{p)+UJr. 



P 



pF^^\p) 
pG^^\p) 



+ m 



pG^^\p) 
-pF^^\p) 



qQ_ 
87r2 



P 



^m{p) 

i^mip) 



Jdq qF(^\q)Q^^. (l^) 
Jdq qG^^){q)Q,^. {^) _ 



qQ_ 
87r2 



q dq 



,<^m(g) 



iF^^'m^i (^) 



qQ_ 
87r2 



m 



qG(^\q)Q^^^ 



\ 2pq J 



4(27r)^ 



dq 



dk 



m 



Q 



P 



kG^^\k)Qj^i 

-fcFW(^)g^.4(^)g^.^.( 



2kq 



2pq 
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kF(^\k)Q^ 



2 \ 2A:q 



2pq 
2pq 



q^Q'^ f dq 
(27r)4 J uJm{q) 



dk < 2m 



A;GW(A;)g, 



2kq 



V'VnJjip^ \^ 2kq J ^JT2 \ 2pq ) 



+q 



where Qj±i is a Legendre function of the second kind and 



(2.15) 



(2.16) 



To rewrite (2.15) in terms of dimensionless variables, the two masses are first 
rewritten as follows: 



m = mo(l-A), M = mo(l + A). 

The dimensionless momentum p' is defined by 



P 



P_ 

mo' 



and the dimensionless energy e by 



E 



E 



e = 



M + m 2mo 



Defining 



co'm(p) 
mo 



= V'(l + A)2+p'^ 



and 



a;_(pO ^ ^ = J(l-A)2+p'2^ 



mo 



(2.17) 



(2.18) 



(2.19) 



(2.20a) 



(2.206) 
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and omitting primes since all variables are now dimensionless, (2.15) becomes 



2eu;+{p) 



pF^^\p) 
pG^^\p) 



;i-A) 



pG^^\p) 
-pF^^\p) 



87r2 



+ 1 



+1^ / g 



[ gGW(,)g,^. (^) 



A) 



'iG^^'m^^, (^) 



{qQ? 



2kq 
2kq 



Q 



2pq 



2pq 



+q 



kF^^\k)Q.., 



\ 2kq 



(27r)4 J u.{q) 



dk{2{l-A) 



kG^^\k)Qj^-^ 
-kF^^\k)Q.^i 



+ q 



kF^^\k)Q^ 



2kq 
k^W 
2kq 



Q 



2pq 



2pq 



VV-«j±i 2kq j -^JTt V 



2pg 



2pg 



(2.21) 
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The set of two equations with the top signs is transformed into the set 
of two equations with the bottom signs and vice versa with the following 
replacements: 

G'(+) ^ F^-\ ^ -G^~\ e ^ -e (2.22) 

As a consequence only the equations for F*^+^and G^^'^ need be solved. For 
notational convenience the superscripts on F^~^^ and G^^^ are omitted in 
future equations. 

3 Solutions to the Nonrelativistic Reduction 
of the Bethe-Salpeter Equation 

In this section the nonrelativistic reduction of the Bethe-Salpeter equa- 
tion, which is just the Schrodinger equation, is solved numerically in mo- 
mentum space. Although the equation can be solved analytically, here it 
is solved numerically to illustrate techniques that will be used to solve the 
Bethe-Salpeter equation in the instantaneous approximation, an equation 
that, apparently, cannot in general be solved analytically. A method is intro- 
duced for handling the singularity in the kernel that makes possible the use 
of basis functions that automatically satisfy the boundary conditions both 
at small and large momenta. For angular momentum states i = and £ = I, 
these basis functions are not an improvement over the basis functions used 
by Spence and Vary [7] that only automatically satisfy the boundary con- 
ditions for small momentum, but for angular momentum states i > 1, the 
basis functions used here converge to a solution much more efficiently. 

Once the instantaneous approximation has been made, it is straightfor- 
ward to make a nonrelativistic reduction [1]. Keeping the lowest order term 
in each of the two interaction terms, proportional to qQ and {qQY, respec- 
tively. 

The nonrelativistic energy E' is related to the relativistic energy E hj E = 
m + M + E' . In the final term in the above equation, the integer "l"in the 
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first parenthesis is from single-plioton excliange wliile tlie integer "-8" is from 
tlie "seagull" interaction. 
Fourier transforming (3.1), 

£;'*(x) = -(— + ^)V'*(x) + — A ^(x) 

+(l-8)(^)2 -i- ^ ^(x). (3.2) 
^ '^An' AM |x|2 ^ ' ^ ^ 

Prom (3.2) it follows that single photon exchange yields a repulsive poten- 
tial term proportional to {qQY that decreases as the square of the distance 
between the constituents while the "seagull" interaction yields an attractive 
potential with the same form that is eight times as strong. To better compare 
the solutions here with those of Spence and Vary [7], in the nonrelativistic 
limit the potential term proportional to {qQY will be neglected. However, 
the effects of this term will be significant when (2.13), the Bethe-Salpeter 
equation in the instantaneous approximation, is solved in the next section. 
To solve (3.1) a dimensionless momentum p' is defined by 

- P (3.3) 



where /i is the reduced mass, /i — Mm/{M + m). Omitting the term pro- 
portional to {qQY, (3.1) becomes 

a+p")*(p')^S^/.:^*W). (3.4, 

Eq. (3.4) is the integral form of the Schrodinger equation for a quanta with 
mass n and charge q interacting with a stationary charge Q via the Coulomb 
potential. Since the momentum variables are all now dimensionless, for no- 
tational convenience the primes will be omitted in future equations. 
The solution is of the form 

vl/(p) = i?(p)F,(^,0). (3.5) 

Using Hecke's theorem [18,19,20,10] the angular integration is easily per- 
formed. The angular dependence separates, yielding an integral equation. 
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- ^ ^r<k QA"^) (3.6) 

Defining 



(3.6) becomes 

{l+p')pR{p) = A- Hdq QeiA^) qR{q). (3.8) 
TT Jo Zpq 

Theoretically, 

X = i + n; n=l,2,.... (3.9) 

Eq. (3.8) would be straightforward to solve numerically were it not for the 
fact that Qeiip"^ + <f)/'^PQ) has a logarithmic singularity at p = g. 

The boundary conditions are determined with the aid of the asymptotic 
relationship for Legendre functions of the second kind [21], 

At small p the function pR{p) is assumed to be of the form 

pR{p) — .p^o (3.11) 

where cq is a constant. From (3.10) it follows that at small p, Qeiip"^ + 
q'^)/2pq) p^^^. Equating the left- and right-hand sides of (3.7), at small p 
the equality pR{p) ~ p^~^^ is obtained, implying 

co = e+l. (3.12) 

At large p the function pR{p) is assumed to be of the form 

pRip)^^. (3.13) 
Using logic analogous to that which lead to (3.12), 

Coo = ^ + 3. (3.14) 
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Solutions are obtained by expanding pR{p) in terms of N cubic B-splines 
BAp) [8], 

N 

pR{p)^F{p)Y,CjBj{p). (3.15) 
i=i 

By choosing the convergence function F{p) in (3.15) so that at small and 
large p it behaves as the solution pR{p) itself, fewer B-sphnes are required to 
represent solutions that go to zero rapidly at the boundaries. 

Cubic B-splines are defined on five consecutive knots. To determine the 
spacing of the knots. A'" — 4 zeros xi of a Chebychev polynomial are calculated 
using the formula 

^i--^s ^2{N^% ^ = 1,2,..., AT -4, (3.16) 

and then the knots 7^+4 on the positive p-axis are determined by 



T,+4 = cJ\^ + C2, i = 1, 2, . . . , AT - 4, (3.17) 

where Ci and C2 are constants. The knot T4 is placed at the origin and 
three knots are placed on the "negative" p-axis to create maximum freedom 
in constructing the solution pR{p) near the origin. The three knots on the 
"negative" p-axis are mirror images of the first three knots in (3.17). 
Spence and Vary [7] note that integrals of the form 



2 2 

dpp'^'QEi^^), fc = 0,l,2,... (3.18) 

2pq 



are both finite and readily calculated analytically, so they choose F[p) in 
(3.15) to be given by F{p) — which, from (3.12), automatically satisfies 
the boundary condition near p = provided that the sum of B-splines in 
(3.15) are non-zero and slowly changing near the origin. Since the B-splines 
vanish at the largest knot, an appropriate sum of B-splines will satisfy the 
boundary condition at large p. However, the boundary conditions are satis- 
fied both at small and at large p with the choice 

(C^ -)-p^j«-|-2 
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where c is a constant. Note that at small p, F{p) p^^^ as expected from 
(3.12), while at large p, F{p) — >• p~^^+'^\ which is one power of p less than is 
indicated in (3.14). Since the last B-spline in the expansion (3.15) vanishes 
at the largest knot, the boundary conditions will be satisfied automatically 
both for small and large momenta provided that the sum of B-splines in (3.15) 
is slowly changing at small momenta and goes to zero as the reciprocal of 
the momentum at large momenta. By choosing F{p) in (3.15) appropriately, 
solutions that decrease rapidly at small and large momenta can be accurately 
reproduced with fewer B-splines than when -F(p) is omitted. 

Eq. (3.8) is solved numerically by converting the integral eigenvalue equa- 
tion into a generalized matrix equation using the Rayleigh-Ritz-Galerkin 
method [21]. The solution is expanded in terms of B-splines using (3.15), 
and then both sides of (3.8) are multiplied by Bi(p)F{p) and integrated over 
p. A generalized matrix equation results that is of the form Ac — XBc where 
the matrices A and B are given, respectively, by 

fOO 

Aij= dpBi{p)F{p){l+p')F{p)B^{p) (3.20a) 

J 

and 

B,j = - / dpB,{p)F{p) / dqQ,{^^^)F{q)B,{q). (3.206) 
TT JQ Jo zpq 

The elements of the column vector c are the expansion coefficients Cj in (3.15). 
Since both of the above matrices are symmetric and Aij is positive definite, 
the eigenvalues arc real [16]. 

The choice F{p) = p^^^ works well for small values of angular momentum 
because the sum of a small number of B-splines readily creates a function 
that decreases asp"^^^"^^^ at large p, thus satisfying the boundary condition as 
given in (3.14). In addition, when F(p) = p^"*"^, the integrals over Quiij? + 
q^)/'2pq) in (3.206) can be performed analytically because they are of the 
form (3.18). For larger values of angular momentum, however, the choice 
F{p) = p^^^ does not work well because the sum of a small number of B- 
splines does not readily create a function that decreases sufficiently rapidly 
at large momentum so as to satisfy the boundary conditions. 

The choice (3.19) for the convergence function immediately allows the 
boundary conditions to be satisfied by a sum of B-splines that is slowly chang- 
ing, but now the integrals over q in (3.206) can no longer be readily performed 



15 



analytically. The method for integrating the variable q over the singularity in 
the integrals is simple conceptually but somewhat involved numerically: Ex- 
cept in an e-neighborhood of the singularity, all integrations are performed 
numerically using Gaussian quadrature with a seven-point option. As the 
integration variable approaches the singularity, where the integrand changes 
most rapidly, integration intervals arc decreased to maintain accuracy of the 
numerical integration. Within an e-neighborhood of the singularity, the in- 
tegrand, excluding the associated Legendre function Qi{{p^ + ?^)/2pg), is 
expanded in a Taylor series about the singularity. The integral is then a 
sum of integrals of the form (3.18) that can be integrated analytically. The 
parameter e is chosen to be the smaller of 0.01 or the distance from the sin- 
gularity to the nearest knot, thus avoiding the complication of integrating 
over a knot. 

To obtain a numerical estimate of the accuracy of each solution, the left- 
and right-hand sides of (3.6) are calculated midway between each pair of 
knots on the (positive) p-axis. A reliability coefficient R [14], which is a 
statistical measure of how closely the two sides of the equation agree at the 
selected points, is calculated. If the two sides of the equation agree exactly 
at all of the selected points, then R equals unity. Determining where the 
left- and right-hand sides of the equation agree least well reveals possible 
problems with trial solutions. 

When F{p) = p^"*"^, excellent solutions are obtained for £ = and i = 
1. With 21 splines in the expansion (3.15), eigenvalues were obtained with 
four or five significant figures, and corresponding i?- values were in the range 
0.999 < R < 1.00. However, when i = 2 as shown in second and third 
columns of Table 3.1, an incorrect eigenvalue appears with a corresponding 
R = 0.00112. For £ = 3, the first few eigenvalues are accurate, but the 
corresponding i?- values are on the order of 10~^, indicating the solutions are 
unreliable. Examination of these solutions reveals that they are incorrect 
at large p. By choosing the convergence factor F{p) as given in (3.18), the 
difficulties that appeared for £ > 1 are eliminated as can be seen from the 
fourth and fifth columns of Table 3.1. 

Table 3.1: Numerical values oi X — £ + 1,£ + 2, ... when 21 splines are used 
in the expansion (3.15). 
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F{p) = p^+i 


F{p) = p^^^ / over{c^ -)-p2j£+3/2 




A 


R 


A 


R 


2 


3.00000 


0.99884 


3.00000 


1.0000 




1.00007 


0.97013 


1.00001 


1.0000 




4.49547 


0.00112 


5.00003 


1.0000 




5.00030 


0.84320 


6.00008 


1.0000 


3 


4.00118 


0.00266 


4.00000 


1.0000 




5.00844 


-0.00028 


5.00000 


1.0000 




6.02762 


-0.00031 


6.00004 


1.0000 




7.06566 


-0.00023 


6.99999 


1.0000 



4 Solutions of the Bethe-Salpeter Equation 
in the Instantaneous Approximation 

Solutions to the Bethe-Salpeter equation in the instantaneous approxi- 
mation arc obtained using two different basis systems. The first basis system 
is comprised of essentially the same basis functions that were employed to 
calculate solutions to the nonrelativistic hydrogen atom. Because the basis 
functions vanish at large momenta, they are particularly suitable for rep- 
resenting solutions that have significant support only to moderately large 
values of momentum (and position). For this basis system four B-splines 
are non-zero between consecutive knots in the physical region except for the 
final four knots at the largest values of momentum: There the number of 
non-zero B-splines between consecutive knots decreases from three to two 
until only one B-spline is non-zero between the final two knots, thus making 
it increasingly difficult to express a solution at large momentum in terms 
of these basis functions. To better represent solutions that are highly local- 
ized in position space and, therefore, have significant support at large values 
of momentum, a second basis system is used in which some basis functions 
vanish only at infinite values of momentum and four B-splines are non-zero 
between consecutive knots in the physical region. 

The boundary conditions as p approaches zero and infinity are determined 
using the same procedure employed in the previous section. The results are 
as follows: 

pG{p) — ^p'^^ pF{p) — >p^^^ (4.1a) 
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Solutions can be obtained using methods of the previous section and are of 
the form 



N 



(4.2a) 



and 



N 



pG{p)=J^,{p)J2fjBj{p). 



(4.26) 



The convergence functions Gi{p) and J^i{p) are chosen so that the boundary 
conditions are automatically satisfied provided that the sums of B-splines in 
the previous equation are slowly changing for small and large momenta. 



Slip) 



pj 



3 + \ 



P 



P 



(4.3) 



In the above equation cq and cjf are constants. Note that at small p, pQi{p) 
and pTiip) vanish as indicated in (4.1a), but at large p they decrease by a 
factor p more slowly than indicated in (4.16) because the B-splines themselves 
vanish at large p. As can be seen from (4.3), solutions go to zero rapidly at 
the boundaries even at the smallest value j = 1/2, so it would be difficult 
to obtain any accurate solutions in the instantaneous approximation without 
using convergence functions. 

Eq. (2.21) is converted into a generalized matrix equation of the form 



9 
f 



= eB 



9 
f 



(4.4) 



by multiplying the top and bottom equations by Gi{p)Bi{p) and J^i{p)Bi{p), 
respectively, and then integrating over p. The elements of the column vectors 
g and / are, respectively, the expansion coefficients gj and fj in (4.2). Since 
the matrices A and B have been constructed so that both are symmetric and 
B is positive definite, the dimensionless energy eigenvalue e is forced to be 
real [16] as required. 
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The Bcthc-Salpeter equation in the instantaneous approximation contains 
double integrals while in the nonrelativistic limit, the equation involves only 
single integrals. In spite of this complication, by performing integrations in a 
specific order, all integrals with a logarithmic singularity that are necessary to 
solve the equation are of the form already encountered in the previous section. 
However, the technique used to integrate over the logarithmic singularity in 
the previous section fails at very large values of momentum: All of the terms 
being expanded in a Taylor series about the singular point p are functions 
of with the result that a typical term in the expansion is ttni^cf — p^)". 
Within an e-neighborhood of p, the maximum value of — is 2pe. When 
p is of the order of 1/e the expansion begins to lose accuracy. By choosing e 
to decrease linearly with increasing p, this problem is avoided. 

To check the accuracy of the solutions by calculating the rehabihty coef- 
ficient, double integrals of the following form must be evaluated: 

r q,^p1+jI) r Q/.l+j^)^tl^B,(k) (4.5) 

Jo u_{q) 2pq ^ Jo ^ 2kq + k^Y^ ' ^ ' 

The integral over the variable k can be calculated as previously discussed. 
Except within an e-neighborhood of the logarithmic singularity of the inte- 
grand, the integral over the variable q is evaluated numerically. Within the 
e-neighborhood, the integral over k is expressed as a power series expansion 
in the variable q, 

foo t-2 _|_ „2\ „£+di 3 

/. * *(^) = - py- (4.6) 

The power series expansion in the above equation depends on the fact that 
the integral vanishes as q^~^^ at small q , a fact that is readily verified using 
(3.10). The coefficients Uj are determined numerically so that the expansion 
and the integral agree atp-|-e,p-|-e/3,p — e/3 and p — e. Using the expansion 
in (4.6), within the e-neighborhood of the logarithmic singularity at p, the 
integral (4.5) can be performed analytically. 

To better represent solutions that are highly localized in position space 
and, therefore, have significant support at large values of momentum, a sec- 
ond basis system is introduced in which some basis functions vanish only at 
infinity. To construct the basis system, the momentum is first mapped onto 
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a compact space with the transformation 

x{p)^b^, (4.7) 

where a and b are constants. Prom (4.7) it follows that —b < x{p) < b. 

The knots are determined by first calculating N — 8 zeros Xi of a Cheby- 
chev polynomial using the formula 

^i--('^^ 2{N-8y ^ = 1,2,..., AT -8. (4.8) 

The knots in the region — 6 < x < 6 are then given by 

Ti+^ = bxi, t = l,2,...,N -8. (4.9) 

The knot T4 is placed at x = —b {p = 0) and three knots arc placed in the 
region x < —b (on the "negative" p-axis) to create maximum freedom in 
constructing solutions near the origin. The three knots in the region x < —b 
are mirror images of the first three knots in (4.9). In a similar fashion the 
knot Tjv-3 is placed at x — b {p — 00) and three knots are placed in the region 
X > b {^^p > 00") to create maximum freedom in constructing solutions at 
very large momenta. The three knots in the region x > b are mirror images of 
the final three knots in (4.9). With the above knot structure, four B-splines 
are non-zero between each pair of adjacent knots in the entire physical region 
< p < 00. 

The solution is expanded in terms of B-splines as follows: 



N 

pG{p) = g2{p{x))J2gjB^{x), (4.10a) 



and 

N 



where 



pG{p) = T2{p{x)) E fjBjix), (4.106) 



p?+2 p)^ 2 

^2(P) = / 2 I 2^7■+l ' •^2(p) = I 2 I 2^7■+2 • (^-11) 
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For the second basis system the final three B-sphnes in the expansion are 
non-zero at x = b {p = oo). Consequently, solutions that have significant 
support at large values of momentum are more readily expressed in terms of 
the second set of basis functions. Since some B-splines are finite at p = oo, 
the functions ^2(p) and J^2{p) ai'e chosen to satisfy the boundary conditions 
(4.1) both at small and large momenta. 

To integrate over singularities at large values of momentum p, the in- 
tegrand is expanded in a Maclaurin series in the variable 1/g instead of in 
a Taylor series. Specifically, if the location of the first knot less than the 
singularity corresponds to a value of momentum equal to or less than 50, 
then integrals are evaluated as previously discussed. On the other hand, if 
the first knot less than the singularity corresponds to a value of momentum 
greater than 50, the integral is evaluated numerically from x — —h {p = Q) 
to the knot. From the knot to x — h {p = cxo), the integral is evaluated 
analytically by expanding the integrand, excluding the Legendre function of 
the second kind, in a Maclaurin series. The necessary formulas for carrying 
out the integration are given in the appendix. 

A corresponding modification is required to evaluate the double integrals 
(4.5). When the location of first knot less than the logarithmic singularity 
at q = p corresponds to a value of p equal to or less than 50, the integral is 
evaluated as before. When the position of the knot corresponds to a value of 
p greater than than 50, the integral is evaluated numerically except within 
an e-neighborhood of the singularity by expanding the integral over A; as a 
Maclaurin series in the variable 1/g, 

/ dk QA T o ,o^^ BAk) = -irrY^H-- (4-12) 

In writing the expansion, the fact has been used that the integral vanishes 
as at large q. The coefficients aj are determined numerically so that 

the expansion and the integral agree at p -|- e, p -|- e/3, p — e/^i and p — e. 
Using the expansion in (4.12), within the e-neighborhood of the logarithmic 
singularity at p, the integral (4.5) can be performed analytically. 

The first basis system has more knots concentrated at small values of 
momentum and, therefore, is more suitable for representing weakly-bound 
solutions or solutions with detailed structure in this region. The second 
basis system has more knots at large values of momentum and is better for 
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representing strongly-bound solutions that have significant support at large 
values of momentum. However, at least when 35 or fewer splines are used, the 
second basis system does not adequately represent the most strongly bound 
solutions with e on the order of or less than about 0.3. Of course, for such 
states the instantaneous approximation is not a satisfactory approximation 
to the Bethc-Salpctcr equation. 

When only the effects of single-photon exchange are included, the solu- 
tions in the instantaneous approximation are approximate solutions of the 
Bethe-Salpeter equation in the ladder approximation. If the masses of the 
constituents are equal, zero-energy, analytical solutions exist for the Bethe- 
Salpeter equation in the ladder approximation for several values of qQ/Aii in 
the range < qQ/An < 100 [10]. In the instantaneous approximation, when 
only single-photon exchange is included, no bound-state solutions were found 
with either zero or finite energy within this range of coupling constants. 

For all data graphed below, the constituent masses are equal although 
solutions with unequal constituent masses are no more difficult to determine. 
Solutions were calculated using the two different basis systems previously 
discussed. For values of e = E/{M -\- m) > 0.95, the graphed results are 
those obtained from the first basis system, which has more knots at small 
momenta. For all other values of e, the graphs are an average of the solutions 
obtained from the two basis systems. Solutions for e obtained from the two 
different basis systems almost always agreed within 0.04 and usually agreed 
more closely while reliability coefficients were almost always greater than 
0.99. Solutions were calculated by expanding the wave function in terms of 
35 B-splines. 

From Fig. 4.1, as the coupling constant decreases in magnitude, the 
repulsive effects of angular momentum become apparent so that states with 
higher angular momentum are more weakly bound. 

Because the bound states in Fig. 4.2 include the "seagull" term, for 
the same value of the coupling constant these states are significantly more 
tightly bound than the corresponding states in Fig. 4.1, which only include 
single-photon exchange. 

The bound states in Fig. 4.3 that occur for positive values of qQ/Ai^ are 
the result of the "seagull" interaction, which is always attractive nonrcla- 
tivistically. As a consequence, like-signed electric charges can attract [10]. 
But, at least for equal- mass constituents in the instantaneous approximation 
as shown in Fig. 4.3(a), no binding occurs between like-signed charges for 
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?g/47r 
(a) 



gg/47r 
(b) 



Figure 4.1: The dimensionless energy eigenvalues e = E/{M + m) of the 
three lowest bound states as a function of the coupling constant qQ/4:7r when 
only single-photon exchange is included in the instantaneous approximation. 
Graphs (a) and (b) correspond to angular momentum j = 1/2 and j = 3/2, 
respectively. 




-0.6-0.5-0.4-0.3-0.2-0.1 -1.2 -1 -0.8-0.6-0.4-0.2 



qQ/^7i qQ/4:7T 

(a) (b) 

Figure 4.2: The dimensionless energy eigenvalues e = E/{M + m) of the 
three lowest bound states as a function of the coupling constant qQ/An in 
the instantaneous approximation. Graphs (a) and (b) correspond to angular 
momentum j — 1/2 and j — 3/2, respectively. 
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gg/47r 
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Figure 4.3: The dimensionless energy eigenvalues e = E/{M + m) of the 
three lowest bound states as a function of the coupling constant qQ/Ar: in 
the instantaneous approximation. Graphs (a) and (b) correspond to angular 
momentum j = 1/2 and j = 3/2, respectively. 



values of qQ/in < 1.9. Since no "elementary" constituent could have such 
a large charge, the effect would not be observable. On the other hand, by 
considering the fully relativistic equation with unequal-mass constituents, 
if such binding could be achieved for constituents that have charges with 
magnitudes on the order of e, the effect might be observable. 

5 Conclusions 

In spite of a derivative coupling in the Lagrangian, the Bethe-Salpeter 
equation, which includes the effects of both single-photon exchange and the 
"seagull" interaction, can readily be solved in the instantaneous approxima- 
tion: Terms appear symmetrically in the equation so that it can be con- 
verted to a generalized matrix eigenvalue equation of the form Ac = eBc 
where the matrices A and B are symmetric and B is positive definite, a suf- 
"Tftient condition for yielding real eigenvalues e. Using a generalization of a 
method introduced by Spence and Vary [7] for handling logarithmic singu- 
larities in integrands, basis systems are used that automatically satisfy the 
boundary conditions, making it possible to obtain solutions in the instanta- 
neous approximation. Weakly-bound solutions are more efficiently obtained 
by writing the solution in terms of basis functions that extend only to finite 
values of momentum while solutions with large binding energies are obtained 
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more accurately by mapping the momentum onto a compact space and then 
expressing the solutions in terms of basis functions, some of which vanish 
only at infinity. A statistical measure is used to provide an indication of the 
accuracy of solutions by determining how well the wave functions satisfy the 
equation midway between each knot in the physical region. 

Acknowledgment s 

I would like to thank Professor John J. Skowronski for very helpful dis- 
cussions regarding statistical parameters that would indicate the reliability 
of numerical solutions to equations. Dr. David G. Robertson assisted in 
optimizing the code. This work was supported by a grant of computer time 
from the Ohio Supercomputer Center. 



Appendix: Calculation of Integrals 

Here the formulas are given for calculating the integrals of the form 
utk, f'^ ^ i V^ + q\ 1 £ = 0; A; = 1,2,3,... , 

that are required to integrate over the logarithmic singularities when they 
occur at large p. Using the recursion formula for Legendre functions of the 
second kind [23], 

a recursion relation for -B^f' ^) (p) follows immediately: 

The integrals -B(f (p) are readily expressed in terms of the integrals 

J a (/ 
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Specifically, 
and 



The integrals I(^a,b)ip) ^{a,b)ip) calculated using standard tables of 
integrals [23], although I(^ab)ip) evaluated as an infinite series. For k > 2, 
the integral I(^ab)^P) calculated using the following formula: Let 

f lri((t + bx) , . , 

I= dx -^—^ k>2. {A7) 

Integrating by parts, 

J. 1 ^ ln{a + bx) , ^ f dx 



k — 1 x''~^ J x^ ^{a + bx) 

The integral in the above expression is evaluated using partial fractions, 
yielding the desired formula: 
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